General approach for studying first-order phase transitions at low temperatures 
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By combining different ideas, a generai and efficient protocoi to deaf with discontinuous phase 
transitions at iow temperatures is proposed. For small T's, it is possible to derive a generic analytic 
expression for appropriate order parameters, whose coefficients are obtained from simple simulations. 
Once in such regimes simulations by standard algorithms are not reliable, an enhanced tempering 
method, the parallel tempering - accurate for small and intermediate system sizes with rather low 
computational cost - is used. Finally, from finite size analysis, one can obtain the thermodynamic 
limit. The procedure is illustrated for four distinct models, demonstrating its power, e.g., to locate 
coexistence lines and the phases density at the coexistence. 
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First-order phase transitions (FOPT) are ubiquitous 
in nature [H, associated to a countless number of pro- 
cesses Q. Moreover, they take place in different temper- 
ature scales, e.g., from 2800 K in the Earth core-mantle 
Q to few and near zero Kelvin range (for many sys- 
tems, in a rather similar way {H), or even being respon- 
sible for unique effects, but across a very broad range of 
T's @. In special, FOPT at low temperatures underpin 
important phenomena, like field-induced metal-insulator 
transitions, magnetoresistance, superfluidity, and Bose- 
Einstein condensation, among many others. 

So, it is understandable and desirable the multitude 
of approaches (mainly numerically @ given the difficulty 
to obtain general exact results |7|) developed to study 
FOPT. In certain instances, nevertheless, many of them 
can face procedural difficulties, not leading to precise re- 
sults for the sough thermodynamical quantities, say, the 
exact location of coexistence lines. 

As particularly powerful methods we can cite cluster 
algorithms [1] , multicanonical @ and Wang-Landau [l(| ■ 
In cluster, non-local configuration exchanges often en- 
sure the crossing of even high free-energy barriers. But 
a drawback is its specialization: each model requires a 
specific and efficient algorithm to implement the transi- 
tions, not available in many cases. On the other hand, 
Wang-Landau and multicanonical are general and have 
been applied successfully in a great diversity of problems. 
However, the former may demand very large computa- 
tional time to calculate the density of states, specially 
considering that the number of states can increases very 
fast with the system size ll|. The latter relies on his- 



togram reweighting techniques to obtain the appropriate 
averages, a difficult task for large systems (see, e.g., [l2|). 



Given so, here we present a protocol to study FOPT 
at low temperatures by means of direct and simple nu- 
merical simulations. Extending previous results |13| . it 
considers around any transition a general parametric an- 
alytical expression for relevant thermodynamic quanti- 
ties (like order parameter, density, magnetization, com- 
pressibility, etc). The parameters are then obtained by 
simulating small systems, making the approach compu- 



tationally fast. Finally, from proper extrapolations, the 
correct thermodynamic limit is obtained. Once standard 
algorithms usually fail for FOPT at low T's, even for 
small systems, we consider tempering methods (already 
shown reliable for FOPT, see 1J,[15[ and Rcfs. therein). 
Thus, we use the parallel tempering, PT, very efficient 
for small and intermediate system sizes. As we exemplify 
with four different lattices models, the approach leads to 
a precise way to determine the coexistence regions. 

We begin recalling a rigorous analysis for finite sys- 
tems having J\f coexisting phases at £*, with £ an ap- 
propriate phase transition parameter control (e.g., tem- 
perature or chemical potential). It has been shown [l6| 
that at low T's and around £*, the partition function 
is very accurately given by Z — J^=i a » e*p[-/3Vf n 
with j3={kT 
free energy [l 
resulting from eventual symmetries of the problem. 

Next observe that W = -d t hx[Z]/{fiV) {d x = d/dx) 
is frequently the start point to calculate distinct order 
parameters (density, magnetization, etc). Since close to 

r, u « /*+/;* v E3 for /;* = d e /„(£*) and y = e-r, 

we find the following general form for W: 



' given uy zj = 2_ 
1 . For the phase n, /„ is the (metastable) 
per volume V and a n is the degeneracy, 



AT 



w = ( & i + 2j bn exp t" 



n2/])/(l + > , °n exp[-a„y]). 



(1) 

The coefficients a n , b n and c„ depend on £*, f'*, T, and 
other system parameters. But only the a„'s are (linear) 
functions of V. Then, at the coexistence (y = 0) W is 
independent on the volume and for all V the curves W x £ 
cross at £ = £*. In this way, Eq. ([T]) can be used not 
only to locate the transition point, but also to determine 
the coexistence order parameters at the thermodynamic 
limit. Moreover, if the /'s are ordered such that /{* = 

J 2 — ■ ■ ■ — Jm ^ Jm+1 J k ~ J fc+1 _ _■ ■ • — J N 

for V -> oo and y ->■ ± we have W± = b n /c n , 
with w + = l,u + = m,V- = fc,M_ = M. For k = M 
(to = I) W+ (W-) is given in terms of the sole phase 
which is immediately to the right (left) of £* . 

Thus, considering relatively small V's we can obtain 
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the parameters a's b's and c's and hence, from the curves 
W, appropriate order parameters and response functions 
- e.g., through derivatives of the order parameter - at the 
transition point. For instance, if £ — p is the chemical 
potential, p(p,T) = —W is the density and \ — <9 m p|t 
is the isothermal compressibility. Finally, from a simple 
scaling analysis 18], but using analytical smooth expres- 
sions, the thermodynamical properties are determined. 

The above will work properly only with methods which 
correctly sample the configuration space [1, [l(| , yielding 
reliable fittings for Eq. ([T]) parameters. Often, this is not 
so when systems displaying strong discontinuous transi- 
tions are simulated by conventional one-flip approaches, 
even for small sizes. The solution is then to consider 
enhanced sampling, like parallel [l9[ and simulated [2(| 
tempering algorithms, PT and ST. Since the former is 
particularly appropriate for FOPT (see 14[ for details as 
well as for implementation), here we use the PT in our 
"combo" procedure for phase transitions at low T"s. 

To illustrate the protocol, next we analyze four differ- 
ent lattice models displaying strong FOPT at low T's. 
In each case, what operationally sets a low temperature 
is the validity of the previous Z decomposition. Physi- 
cally, it corresponds to T's for which there is no overlap 
between the peaks of the order parameter bimodal dis- 
tribution at the coexistence. In all examples we perform 
accurate numerics with the PT algorithm and compare 
with the general Eq. ([TJ, whose parameters are always 
obtained using only four points from the simulations. 

As the first example, we consider a rather complex sys- 
tem, the associative lattice-gas (ALG) model [23|, aimed 
to reproduce liquid polimorphism and water-like anoma- 
lies. A site i may or may not be occupied [a t — 1 or 0) 
by a molecule in a triangular lattice. The orientational 
variable T t J — 0, 1 represents the possibility of hydrogen 
bonding (in a maximum of four) between the molecule in 
site i and those in the adjacent six neighbors j, provided 
T i J T j l = !• Two first neighbor molecules have an inter- 
action energy of — v (— v + 2u) if there is (there is not) a 
hydrogen bond between them. The Hamiltonian reads 



U = 2u 



i [(l-v/(2u)) 



T T 

1 3 



<i,3> 



It presents one gas and two liquid, LDL and HDL, phases 
of densities p = 0, p = 3/4, and p = 1, respectively. For 
fixed T, by increasing p we pass through two FOPT, 
namely, gas-LDL and LDL-HDL. 

We study the ALG model gas-LDL and LDL-HDL 
FOPT by plotting the density vs p for T = 0.300, 
u = v = 1, and different Vs. For the gas-LDL case, 
we show the results in Fig. Q] (a) . We clearly see a good 
coincident crossing of all curves at p = —1.9986(2), for 
p = 0.600(1). The exact transition density p = 3/5 = 0.6 
is understood recalling that at the coexistence both gas 
(p = 0) and LDL (p = 3/4) phases have equal weight. 
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FIG. 1. The ALG model for the parameters as in the text and 
different V — L x L. The continuous lines are curves from 
Eq. {1} for: (a) p X p around the gas-LDL transition (the 
inset shows p x /i in the whole /i range where the two FOPS 
take place); and (b) </)(= 4p — 3) x fj, around the LDL-HDL 
transition. The blow up shows the crossing for L = 8, 12, 20. 
The expected linear a x V behavior is displayed in the inset of 
(b) (the log-log scale is just for accommodating both cases). 
For comparison, the simulations for V = 4 x 2 are exact. 



Given that a^DL = 4, the value follows. Around LDL- 
HDL, p does not vanish, inset of Fig. [T](a). Since 3/4 (the 
totality) of the lattice is filled by molecules in the LDL 
(HDL) phase, a better order parameter is the rescaled 
density 4> = (4,0 — 3). Thus, in Fig. [2(b) we display 4>x p 
for the LDL-HDL transition. Again, all the isotherms are 
well described by Eq. (JT]), crossing at p = 2.0000(3) with 
p = 0.857(1). In Fig. Q](b) inset we confirm the expected 
linear dependence on V for the parameter 02 = a (a sin- 
gle a once we have only two phases in each transition). 

Next we consider the Bell-Lavis (BL) model, which also 
displays water-like anomalies. The sites may or may not 
be occupied (a t — 1 or 0) by molecules of two possible 
orientations. But differently from the ALG, the van-der- 
Waals interaction between two adjacent molecules is at- 
tractive, —e v dw So, there is no energetic punishment if 
hydrogen bonds (of energy —thb) are not formed. Such 
distinctions from the ALG, e.g., result in a second order 
phase transition for LDL-HDL, but still a FOPS for the 
gas-LDL. It is described by (C = e vdw /£hb) 



n 



ai a i ^ 



v 31 
3 



•C] 



(3) 



<»j> 



For C < 1/3, the BL presents three phases, gas (p = 0), 
LDL (p = 2/3 in a honeycomb structure), and HDL (p = 
1) [21(. In the numerics we set £hb = 1 and e v d w = 1/10. 

For T = 0.25, in Fig. [2]we plot pxp around the FOPT 
gas-LDL. Once more, the simulations are well described 
by Eq. (JXJ) . The isotherms cross at p = —1.6528(1), with 
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FIG. 2. For the BL model gas-LDL FOPT and parameters 
as in the text, px/i for distinct Vs and T = 0.25. All the 
continuous lines are properly obtained from Eq. {TJ. The 
upper insets show the isothermal compressibility \ vs p and 
pv (the values of p at the peaks of x) vs V~ . The lower inset 
shows p xT curves for p — —1.6528, which cross at T = 0.25. 

p « 0.507(2) very close to the exact p = 1/2 (which can 
be inferred as done for the ALG model). In the upper- 
left inset we show \ — (§^)t by properly differentiat- 
ing Eq. ([1]) (continuous lines) and by numerically sim- 
ulating x = V((p 2 ) ~~ (p) 2 )- Note the remarkable agree- 
ment, again illustrating the power of Eq. ((TJ) to describe 
relevant thermodynamic quantities around FOPT. The 
upper-right inset displays the values of p = pv (for which 
X is maximum) vs V . This type of scaling extrapola- 
tion also can give the thermodynamic limit for the transi- 
tion, here p = —1.6527, basically the same value obtained 
from the crossing. Finally, instead of p one could take 
T as the control parameter. Setting p = —1.6528 and 
varying T we see in the lower inset of Fig. [5] the gas-LDL 
transition. As it should be, the curves cross at T = 0.25, 
with p « 1/2. Finally, we note that for T > 0.43, the re- 
sults from the present method starts to be less accurate. 
The Blume-Emery-Griffiths (BEG) model yields HI 

n = - E i Ja * °i + K °i°i] - E i Ha * - Da ^ ( 4 ) 

<i,j> i 

where a site i is either empty or occupied by two dif- 
ferent type of species (cjj = 0,±1). Parameters J and 
K are interaction energies and D and H denote linear 
combination of the species /x's. This system is a particu- 
larly interesting test because the otherwise very reliable 
cluster algorithm for the BEG model § fails for some 
particular K/J's, e.g., the value we address K/J = —0.5. 
So, for a better comparison with our procedure, we also 
propose a new cluster-Metropolis hybrid approach, which 
includes intermediary Metropolis algorithm steps (details 
will appear elsewhere). We note, nevertheless, that the 
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FIG. 3. For the BEG model, p x p for parameters as in the 
text and distinct Vs. The results for V = 4 x 4 are exact. 
The continuous lines properly come from Eq. JTJ. The right 
inset shows Dv {D for which x is maximum) vs V 1 . The 
left inset compares numerical simulations from the PT, cluster 
and hybrid cluster-Metropolis algorithms for V = 10 x 10. 

Metropolis alone is not able to cross the high free energy 
barriers at the phases coexistence. 

In Fig. M we plot p x D for H = 0, J = 1, K = -0.5, 
T = 0.20, and different Vs. We have a FOPT with all 
the isotherms crossing at D — 0.9984(1) and p « 2/3. 
The right inset of Fig. [3] shows the position Dy of the 
peak of x — ~(§d^ t ' calculated directly from Eq. ([1]). A 
linear extrapolation of Dy x V~ x gives D — 0.99845(5), 
in excellent agreement with the crossing value. For 
l/=10xl0, we plot in the left inset simulations from 
the usual cluster, the improved (but dedicated) cluster- 
Metropolis, and PT algorithms. The latter two display 
very good concordance, with the cluster given poorer re- 
sults. We should mention that for the BEG and ALG 
models there are no precise simulations in the literature 
for the parameter conditions here considered. 

Lastly, we discuss the asymmetric Ising Hamiltonian 
on a triangular lattice (of sublattices (a,/3,7)) [25[ 

ti = — j <Ti<jj - k ^2 vt^jCk - .gyVy (5) 

<i,j> <i,j,k> i 

The second sum is over first neighbors trios forming tri- 
angles. Using the Wang-Landau method the model 
has been studied in details (2|| (but for larger systems 
and lower numerical precision). It displays one ferri- 

magnetic, ( h), and two ferromagnetic, (+ + +) and 

( ), phases. For very low temperatures, by increas- 
ing H the system displays a second-order phase transition 

( ) -> (--+) and then a FOPT (--+) ->• (+ ++)• 

The (— — +) phase disappears in a critical endpoint 
(T c = 2.443(1), H c = -2.934(1)), above it giving rise 
only to a FOPT between the two ferromagnetic phases. 
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FIG. 4. For the asymmetric Ising model and parameters as 
in the text, m x H for distinct Vs and T = 5.00 (T = 2.40, 
left inset). Continuous lines are the curves from Eq. |TJ. The 
right inset shows P m x m for the T — 5.00 case. 



Although the magnetization per site m is not the actual 
order parameter, for rather small system sizes we can 
extract from it any relevant FOPT information. 

In Fig. H] we plot m X H for J = 1, K = 2 and 

T = 5.00 for the ( )->(+ + +) FOPT. We see 

that Eq. (JTJ represents quite well the transition. In the 
right inset we show the histogram magnetization den- 
sity probability P m vs m for L = 12, T = 5.00 and 
ff = —2.3325, illustrating further that the phases coexis- 
tence is being properly characterized 14}. Likewise, the 
FOPT ( h) -)■ (+ + +) for T = 2.40 in the left inset 



is well described by our method. All the isotherms cross 
at H = -2.2863(5) and H = -2.9357(5) (left inset), in 
fair agreement (given the different numerical accuracies) 
with the estimates H = -2.284(1) and H = -2.939(1) 
by the authors of Ref . [26j (private communication) . 

By considering Eq. ([IJ, derived from rigorous results 
at low T's, we have proposed a general protocol to study 
FOPT. It is accurate and demands only few simulations 
for relatively small systems, hence a computationally low 
cost procedure. The approach has been very successfully 
applied to four distinct lattice models. Of course, more 
analyzes, e.g., for higher dimensions and continuous sys- 
tems (presently under progress, but with promising pre- 
liminary findings) are in order as further tests. Neverthe- 
less, we believe the method already shows itself a valuable 
tool to analyze the very important problem of FOPT at 
low temperatures. 
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